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Abstract 

We expose in full detail a constructive procedure to invert the so-called "finite 
Markov moment problem" . The proofs rely on the general theory of Toeplitz ma- 
trices together with the classical Newton's relations. 
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Resume 



Nous presentons en detail une procedure constructive pour inverser le 
"probleme fini des moments de Markov". Les preuves reposent sur la 
theorie generale des matrices de Toeplitz et les classiques relations de 
Newton. 
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Version francaise abregee 

Afin d'inverser le systeme fini et mal conditionne (1), Koborov, Sklyar 
et Fardigola, [9,13] ont propose un algorithme recursif non-lineaire. Dans 
[7] nous avons prouve un Icmmc Ic rcduisant a I'extraction de valeurs 
propres generalisees, voir (4). Ccttc Note vise a expliquer en detail les 
raisons pour lesquelles cette procedure simplifiee resout le probleme des 
moments de Markov. Apres avoir rappelc quelques elements de la theorie 
des matrices de Toeplitz et les relations de Newton (Propositions 1 et 
2), nous reformulons cet algorithme simplifie afin d'etablir facilement 
certains lemmes techniques. Finalement, le Thcoreme 1 demontre le lien 
entre valeurs propres generalisees (4) et I'inversion de (1). 



1 Introduction 

We aim at inverting a moment system often associated with the presti- 
gious name of Markov, as appearing in [2,3,6,11,12] in several fields of 
application; consult [5,10,14] for general background on moment prob- 
lems. The original problem is the following. Given the moments ruk for 
k = 1, . . . ,K, find a bounded measurable density function / and a real 
value X > such that 

(1) lo"" fiOe-'d^ = mk, k = l,...,K, 

(2) 1/(01 = 1 almost everywhere on ]0,X[, 

(3) / has no more than K — 1 discontinuity points inside ]0, X[. 

The solution is a piecewise constant function taking values in {—1, 1} a.e. 

on ]0,X[ and changing sign in at most K — 1 points, which we denote 
{uk}, ordered such that < -Ui < ■ ■ ■ < = X. Finding {u^} from 
{ruk} is an ill-conditioned problem when the Uk values come close to 
each other; its Jacobian is a Vandermonde matrix and iterative numeri- 
cal resolution routines require extremely good starting guesses. For less 
than 4 moments, however, a direct method based on solving polynomial 
equations was presented in [12]. Here we are concerned with an arbitrary 
number of moments K eN. 

We consider a slightly modified version of the problem which is more 
relevant for us. In this case / takes values in {1,0} instead of { — 1,1} 
and the moments are scaled as kirik. Moreover, to simplify the 

discussion we confine ourselves to the case when K is even, setting K =: 
2n. The resulting problem can then be written as an algebraic system of 
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nonlinear equations: Given ruk find Uk such that 



rrik 



E 



u 



2j 



U 



(1) 



An algorithm for solving this problem was presented by Koborov, Sklyar 
and Fardigola in [9,13]. It requires solving a sequence of high degree 
polynomial equations, constructed through a rather complicated process 
with unclear stability properties. In [7] we showed that the algorithm 
can be put in a more simple form that makes it much more suitable for 
numerical implementation. The simplified algorithm reads 

(1) Construct the matrices A and B: 



A 



-nil 



-m2„-i . . . -mi 2n 



B 



mi 



^^m2„-i ... mi 2n^ 



(2) 



(2) Let m — (mi, m2, . . . , m2nY solve the lower triangular Toeplitz 
linear systems 

Aa — m, Bb — —m, (3) 

to get a and h. 

(3) Construct the matrices Ai, A2 from a — (oi, 02, ... , a2n)'^ as 



Ai 



( \ 

01 02 ... a„ 

02 03 ... Cln+l 



02n-l 



/ \ 

02 03 ... a„_|_i 

03 04 ... a„+2 



^a„+i a„_|_2 ■ ■ ■ 02n y 



and the corresponding matrices i?i, i?2 from h. 
(4) The values {wfe} can then be computed as the generalized eigenvalues 
of the problems 



A2V — uAxV^ B2V — uBiV. 



(4) 



The forthcoming section is devoted to a complete justification of this 
algorithm. We recall that these inversion routines have been shown to be 
numerically efficient in the paper [7]. 



3 



2 Analysis of the algorithm 



We begin by stating two classical results of prime importance for the 
analysis. 

Let L" C M"^" be the set of lower triangular n x n real Toeplitz matrices. 
We define the diagonal scaling matrix and the mapping T : M" — > L" as 







n-1 



( 



T{x) 



Xi 
X2 Xi 



\ 



^x^ 

X2 



X 



y 10 — ± J \^Xn ... X2 X\j yx„ y 

The mapping T has the following properties, see e.g. [1]: 

Proposition 1 Lower triangular Toeplitz matrices commute and L" is 
closed under matrix multiplication and (when the inverse exists) inver- 
sion, 

T{x)T{y) = T{y)r{x) e L", T{x)-' e L" (5) 

Moreover, T is linear and 



T{x)y = T{y)x, T{x)T{y) = T{T{x)y). 



(6) 



The A matrix has the property 



T{Ax) = AT{x) - T{x)A. 



(7) 



Another result that we rely heavily on is the classical Newton relations, 
see e.g. [8]: 

Proposition 2 (Newton's relations) Let P be the n-degree polyno- 
mial, 

P{x) = Co + CiX H h CnX"' =: Cn{x — Xi) ■ ■ ■ {x - Xn). 

Set Sq = n and define Sk for k > as the sum of the roots of P taken to 
the power k, Sk — S"=i a;^ . Then, the n-\- 1 following relations hold: 

CkSo + Ck+iSi H h CnSn-k = kck, /c = 0, . . . , n. (8) 
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2. 1 Reformulation of the simplified algorithm 



We want to write the equation Aa — m using the mapping T: hence we 
augment the m and a- vectors with a zero and one clement, respectively, 
to get rfi = (0, m)-^ and a = (1, a)'^, both in M-^^-'^. We observe that the 
74-matrix in (2) is the lower right K x K block of A — T{rn). Therefore, 



(A -T{rh))d = -rh + 




( 

V 







^Aa — m 

Thus the equation Aa — mm (3) is equivalent to 

T{rh)a = Aa. 

By the same argument, Bb — —m in (3) is equivalent to 

T{rh)h = -A6. 

with b = (1, bf. 



(9) 



(10) 



We can then directly also show that T(d) and T(6) are in fact each 
other's inverses. 

Lemma 1 When a are b given by (9, 10) then T{a)T{b) — I. 
Proof: By Proposition 1 and (9, 10), 



AT(a)6 = T(Aa)6+T(a)A6 = T(T(m)a)6+T(a)A6 = T(a) 



T(rh)b + Ab 



Thus T{a)b lies in the nuUspace of A which is spanned by the vector 
1 := (1. 0, 0)"^. Moreover, {T{a)b)i = 1 since the first elements of a 
and b are both one and we must in fact have T(a)6 = 1. The lemma 
then follows from (6) and the definition of T. □ 



2.2 What lies beneath the algorithm 



To understand the workings of the algorithm we need to introduce some 
new quantities and determine how they relate to a, b and rh. 

Let us start with some notation: we set Xj — U2j-i and yj — U2j for 
j — 1, ...,n. Furthermore, we introduce the sums 



n n 

^fe = E4' ^*^ = E?/^ k = l,2,...,K = 2n, (11) 



and define Xq — Yq — K. In the even case, it then holds that 



mfe = - = Xfe - Ffc, A; = 1, . . . , X = 2n. 

i=i i=i 

We also define the two polynomials 

p{x) = {x - Xi) ■ ■ ■ {x - Xn) =: Co + CiX -\ h Cn-lX"'~^ + CnX"^ , (12) 

and 

q{x) ^ {x - yi) ■ ■ -{x - Un) -.dQ + d^x^ h dn-ix"^'^ + dnx"". (13) 

We note here that by construction c„ = (/„ = 1. 
By applying (8) to x^p{x) with A; = 0, . . . , K we get 



^0 

Xi Xq 



Cn-1 



^Xk . . . Xi Xq j 



Co 






V ^ y V 



KCn 

{K - l)c„_i 

{K - n)co 






The analogous system of equations holds also for and d^. We introduce 
now some shorthand notation to write these equations in a concise form. 
First we set c = (cn, . . . , cq)^ G M"+^ and d = (o?„, . . . , o?o)^ e ]R"+^ 
We then construct the larger vectors, padded with zeros: c = (c, 0)^ 
and d = (d, 0)'^, both in R^+^. Finally we let X = (Xq, . . . , XkV and 
1^ = (Yq? • • • ) ^k)^ ■ Using T and A we can state the systems of equations 
above as follows: 

T{X)c = {KI - A)c, T{Y)d = {KI - A)d. (14) 

We also clearly have rh = X — Y . 

Before we can relate c and d with d we need the following lemma: 

Lemma 2 Let / : ^ M" &e defined by f{x) := T{x)^^Ax for x with 
a non-zero first element. Then f{xi) = f{x2) implies that Xi = ax2 for 
some non-zero ck e R. 
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Proof: Suppose f{x) = y. Then T{x)y = Ax and by Proposition 1, 
{T{y) — A)x = 0. Hence f{xi) = f{x2) implies that Xi and X2 both 
he in the nuUspace of T{y) — A. Since the top left element of A is zero, 
it follows that the first element of y is zero and therefore the diagonal 
of T{y) is zero. Consequently, the nuUspace of T{y) — A has the same 
dimension as that of A, which is one. □ 

We can now merge together and express the general structure from (12, 
13) and (9, 10) in the most concise way. 

Lemma 3 Suppose c, d are defined by (12, 13) and d, b by (9, 10). 

Then 

T{d)c = d, T{h)d = c. 

Proof: We only need to prove the left equality. The right one follows 
immediately from Lemma 1. Let v = T(d)c = T(c)d. We want to show 
that V = d. We note first that by (14) 

T{X)c^{KI-A)c =^ T{c)X ^{KI-A)c X ^ KT{c)-^c-T{c)-^Ac, 

where T(c) is invertible since c„ = 1. Moreover, it is clear that T{y)l — y 
for all y. Hence, X — Kl — T(c)~^Ac. In the same way we also obtain 
Y^Kl- T{d)-^Ad. Then, 

r(m)d = r(d)m = T{d){X - Y) 

= -r{d)T{c)-^Ac + T{d)T{d)-^Ad. 

We now note that by Proposition 1, 

T{d)r{c)-^Ac = T{c)-^T{d)Ac = T{cy^T{Ac)d = T{c)-^AT{c)d - Ad 
= T{c)-^Av - Ad. 

Since also, T{c)T{d) = T{v) we get 

T(m)d = -T{c)-^Av + Ad + T{d)T{d)-^Ad 



r(d) 



T{d)-^Ad-T{v)-^Av 



+ Ad. 



Consequently, by (9), T{d)~^Ad = T{v)~^Av and by Lemma 2, v = ad] 
for some a e M. But for the first element in v we then have vi = c„ = adn 
and we get a — 1 since Cn — dn — 1. □ 

Finally, we also establish the following lemma. 

Lemma 4 Let V and W be the Vandermonde matrices corresponding to 
the roots of p{x) :— xp{x) and q{x) :— xq{x) respectively. Then 

V^RT{c)V = diag({p'(a;fc)}), W^RT{d)W = dmg{{q'{xk)}), 
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where R — {5n+2-i-j} £ ^5 the reversion matrix and xq — yo — 

0. 

Proof: We have 

n r 

{V^Rr{c)V),, = xt,RT{c)x,., = J2J2crxUxy_i = 

r=0 e=o 

showing the left equahty. The right equahty follows in the same way. □ 



p{Xi-l)-p{Xj-l) iJ^j 



2.3 Conclusion 



We can now conclude and show that the unknown values Uj in (1) are 
indeed the generalized eigenvalues of (4). 

Theorem 1 Suppose K — 2n; let a, b be defined by (3). If all values 
{xj} U {yj} are distinct, then {xj}, {yj} are the generalized eigenvalues 

of (4). 

Proof: Let a, b be defined by (9, 10), which is equivalent to (3). Also 
define c and d as before by (12, 13). By Lemma 3 we have T{a)c -de 









Oi 1 






a„ ... Oi 


1 






ai 


1 


O'n+2 . . . 




Oi 1 


'^271 02n-l • • • 




ttn-l . . . ai 1 
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d 





\ / 



(15) 



Clearly, the lower left block of the matrix multiplied by c is zero, i.e. 
Z^iLo'^i^i+fc = for /c = l,...,n. Now, let Vi be the coefficients of the 
polynomial v{x) :— p{x)/{x — Xj) for some fixed j. Hence, by the special 
structure of (12), 



Co + C\X + ... + CnX^ —. {Vi + V2X + ... + VnX^ ){x - Xj), 
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and, for i — 0, n, 



Ci = < 



-XjVi+i, i = 0, 

Vi — XjVi+i, 1 < i < n — 1, 
Vi, i — n. 



Thus we deduce, 



n— 1 



i=l 



i=l 



i=l 



which is the componentwise statement of A2V = XjAiV. It remains to 
show that the rightmost sum is non-zero for at least some k, so that Xj 
is indeed a well-defined generalized eigenvalue. Let a = (1, 
and V — {0,vi, . . . , . Then (15) gives T{a)c — d and, using Lemma 
4 while taking A; = 1 we have the sum 



i=l 



since V^Rd = {q{xk)} and V^v = {xkv{xk)} = {Sk-jP'{xk)}. Hence, the 
sum is non-zero because Xj 7^ yi for all i,j. The same argument can be 
used for any j, which proves the theorem for {xj}. The proof for {i/j} is 
identical upon exchanging the roles of c, d and d, b. This leads to (4). 
□ 
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